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Abstract We investigate the dynamical evolution of 
hierarchical three-body systems under the effect of tides, 
when the ratio of the orbital semi-major axes is small 
and the mutual inclination is relatively large (greater 
than 20°). Using the quadrupolar non-restricted ap- 
proximation for the gravitational interactions and the 
viscous linear model for tides, we derive the averaged 
equations of motion in a vectorial formalism which is 
suitable to model the long-term evolution of a large 
variety of exoplanetary systems in very eccentric and 
inclined orbits. In particular, it can be used to derive 
constraints for stellar spin-orbit misalignment, capture 
in Cassini states, tidal-Kozai migration, or damping of 
the mutual inclination. Because our model is valid for 
the non-restricted problem, it can be used to study 
systems of identical mass or for the outer restricted 
problem, such as the evolution of a planet around a 
binary of stars. Here, we apply our model to three dis- 
tinct situations: 1) the HD 80606 planetary system, for 
which we obtain the probability density function dis- 
tribution for the misalignment angle, with two pro- 
nounced peaks of higher probability around 53° and 
109°; 2) the HD 98800 binary system, for which we 
show that initial prograde orbits inside the observed 
disc may become retrograde and vice-versa, only be- 
cause of tidal migration within the binary stars; 3) the 
HD 11964 planetary system, for which we show that 
tidal dissipation combined with gravitational perturba- 
tions may lead to a decrease in the mutual inclination, 
and a fast circularization of the inner orbit. 
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1 Introduction 

At present, about 50 multi-planet systems have been re- 
ported, out of which roughly 1/3 possess close-in plan- 
ets with semi- major axis smaller than 0.1 AU. There are 
also indications that about half of these stars are likely 
to have distant companions (Fischer et al, 2001). In 
addition, close binary star systems (separation smaller 
than 0.1 AU) are often accompanied by a third star (e.g. 
D'Angelo et al, 2006). Therefore, although hierarchical 
systems are considerably different from our Solar sys- 
tem, they represent a significant fraction of the already 
known systems of stars and planets. These systems are 
particularly interesting from a dynamical point of view, 
as they can be stable for very eccentric and inclined or- 
bits and thus present uncommon behaviors. In partic- 
ular, they become interesting when the two innermost 
bodies are sufficiently close to undergo significant tidal 
interactions over the age of the system, since the final 
outcome of the evolution can be in a configuration that 
is totally different from the initial one. 

The origin and evolution of the orbital configura- 
tions of multi-body systems can be analyzed with direct 
numerical integrations of the full equations of motion, 
but the understanding of the dynamics often benefits 
from analytical approximations. Additionally, tidal ef- 
fects usually act over very long time-scales and there- 
fore approximate theories also allow to speed-up the nu- 
merical simulations and to explore the parameter space 
much more rapidly. Secular perturbation theories based 
on series expansions have been used for hierarchical 
triple systems. For low values of the eccentricity, the ex- 
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pansion of the perturbation in series of the eccentricity 
is very efficient (e.g. Wu and Goldreich, 2002), but this 
method is not appropriate for orbits that become very 
eccentric. An expansion in the ratio of the semi-major 
axis a\jai is then preferred, as in this case exact ex- 
pressions can be computed for the secular system (e.g. 
Laskar and Boue, 2010). 

The development to the second order in a\/a,2, called 
the quadrupole approximation, was used by Lidov (1961, 
1962) and Kozai (1962) for the restricted inner prob- 
lem (the outer orbit is unperturbed). In this case, the 
conservation of the normal component of the angular 
momentum enables the inner orbit to periodically ex- 
change its eccentricity with inclination (the so-called 
Lidov-Kozai mechanism). There is, however, another 
limit case to the massive problem, which is the outer re- 
stricted problem (the inner orbit is unperturbed) . Palacian 
et al (2006) have studied this case and discussed the ex- 
istence and stability of equilibria in the non-averaged 
system. Farago and Laskar (2010) derived a simple model 
of the outer restricted case and described the possible 
motions of the bodies. They also looked at the quadrupo- 
lar problem of three masses and show how the inner and 
outer restricted cases are related to the general case. 

For planar problems, the series expansions in a\jai 
should be conducted to the octopole order (e.g. Mar- 
chal, 1990; Ford et al, 2000; Laskar and Boue, 2010), 
as the quadrupole approximation fails to reproduce the 
eccentricity oscillations (e.g. Lee and Peale, 2003). How- 
ever, the inclinations of the already known hierarchical 
systems have not been yet determined, and it can be as- 
sumed that high values for the eccentricities may also 
indicate that their mutual inclinations are large as well 
(e.g. Laskar, 1997; Chatterjee et al, 2008). 

As for Mercury, Venus and the majority of the nat- 
ural satellites in the Solar system, close-in bodies un- 
dergo significant tidal interactions, resulting that their 
spins and orbits are slowly modified. The ultimate stage 
for tidal evolution is the synchronization of the spin 
and the circularization of the orbit. Indeed, the ob- 
served mean eccentricity for planets and binary stars 
with a\ < 0.1 AU is close to zero within the obser- 
vational limitations (e.g. Pont et al, 2011). Although 
tidal effects modify the spin in a much shorter time- 
scale than they modify the orbit, synchronous rotation 
can only occur when the eccentricity is very close to 
zero: the rotation rate tends to be locked with the or- 
bital speed at the pcriapsis, because tidal effects arc 
stronger when the two bodies are closer to each other. 

During the formation process, the orbital eccentric- 
ity can increase due to gravitational scattering, so that 
the inner bodies become close enough at periapsis for 
tidal interactions to occur (e.g. Nagasawa et al, 2008). 



The same gravitational scattering is simultaneously re- 
sponsible for an increase of the mutual inclination of 
the orbits (e.g. Chatterjee et al, 2008), and the fact 
that inclined systems exchange its inclination with the 
inner's orbit eccentricity, results that the dissipation 
in the eccentricity can be transmitted to the inclina- 
tion of the orbits, and vice-versa. The most striking 
example is that the spin and the orbit can be com- 
pletely misaligned (e.g. Pont et al, 2009; Triaud et al, 
2010). Previous studies on this subject have been un- 
dertaken by Eggleton and Kiscleva-Eggleton (2001) for 
binary stars and by Wu and Murray (2003) and Fab- 
rycky and Tremaine (2007) for a planet in a wide bi- 
nary. Despite the success obtained by these works in ex- 
plaining the observations, they all used the same set of 
equations, derived by Eggleton and Kiseleva-Eggleton 
(2001), which is not easy to implement and has been ob- 
tained in the frame of the inner restricted quadrupolar 
approximation. As a consequence, their model cannot 
be applied to a large number of situations, where the 
outer orbit cannot be held constant, such as regular 
planetary systems or planets around close binaries. 

In this paper we intend to go deeper into the study 
of hierarchical three-body systems, where the inner- 
most bodies undergo tidal interactions. We do not make 
any restrictions on the masses of these bodies, and use 
the quadrupolar approximation for gravitational inter- 
actions with general relativity corrections. Our study 
is then suitable for binary star systems, planetary sys- 
tems, and also for planet-satellite systems. We also 
consider in our model the full effect on the spins of 
the two closest bodies, including the rotational flatten- 
ing of their figures. This allows us to correctly describe 
the precession of the spin axis and subsequent capture 
in Cassini states. We adopt a viscous linear model for 
tides (Singer, 1968; Mignard, 1979), as it provides sim- 
ple expressions for the tidal torques for any eccentricity 
value. Since we are interested in the secular behavior, 
we average the motion equations over the mean anoma- 
lies of the orbits and express them using the vectorial 
methods developed by Boue and Laskar (2006), Correia 
(2009), and Tremaine et al (2009). 

In Section 2 we derive the averaged equations of mo- 
tion that we use to evolve hierarchical systems by tidal 
effect. In Section 3 we obtain the secular evolution of the 
spin and orbital quantities in terms of reference angles 
and elliptical elements, that are useful and more intu- 
itive to understand the outcomes of the numerical sim- 
ulations. In Section 4 we apply our model to three dis- 
tinct situations of extra-solar systems: HD 80606, HD 98800, 
and HD 11964. Finally, last section is devoted to the 
conclusions. 
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2 The model 

We consider here a hierarchical system of bodies com- 
posed of a central pair with masses mo and to-i , together 
with an external companion with mass TO2- Both inner 
bodies are considered oblate ellipsoids with gravity field 
coefficients given by Ji a and J^ , rotating about the axis 
of maximal inertia along the directions §o and §1 (gyro- 
scopic approximation), with rotation rates loq and u>i, 
respectively, such that (e.g. Lambeck, 1988) 



■h; = h 



3Gto, 



(1) 



where G is the gravitational constant, Ri is the radius 
of each body, and k2 i is the second Love number for 
potential (pertaining to a perfectly fluid body). 

We use Jacobi canonical coordinates, with ri being 
the position of toi relative to too, and r? the position 
of TO2 relative to the center of mass of too and mi . We 
further assume that |ri| <IC |r*2 1 , and we shall refer to 
the orbit of to-i relative to too as the inner orbit, and 
the orbit of to-2 relative to the center of mass of too and 
mi as the outer orbit. In the following, for any vector 
u, u = u/ ||u|| is the unit vector. 



2.1 Conservative motion 

In the quadrupolar three-body problem approximation, 
the potential energy U of the system is given by (e.g. 
Smart, 1953): 
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where Pzix) = (3x 2 — l)/2 is the Legendre polynomial 
of degree two, and terms in (ri/r2) 3 and (Ri/rj) 3 have 
been neglected (i, j = 0, 1). We also have m i = (to + 
toi), fii = to toi/to i, f3 2 = to ito 2 /(to i + to 2 ), /«i = 
Gto 01 , and fi 2 = G(m 01 + to 2 ). 

The evolution of the spins can be tracked by the 
rotational angular momenta, Lj ~ CiUJi §,. In turn, the 
evolution of the orbits can be tracked by the orbital an- 
gular momenta, Gj = fii^J HiCLill — ef) kj (where kj is 
the unit vector G,), and the Laplace- Runge-Lenz vec- 



tor, which points along the major axis in the direction 
of pcriapsis with magnitude e\\ 

*i xG i r i ( ^ 

ei = — s ■ (3) 

Pi mi n 

cii is the semi-major axis (that can also be expressed 
using the mean motion, m = y/ni/a 3 ), et is the eccen- 
tricity, and Ci is the principal moment of inertia. The 
contributions to the orbits are easily computed from 
the above potentials as 



Gi = ri x Fi , G 2 — r 2 x F 2 , 

and 

1 ( G i • ^ 

ei = Fi x — - + f i x Gi I 

piMi V pi J 



(4) 



(5) 



where F 4 = -V Ti U', with V = U + Gm mi/ri + 
GTO iTO 2 /r 2 . 

In Jacobi coordinates, the total orbital angular mo- 
mentum is equal to Gi + G2 (e.g. Smart, 1953). Since 
the total angular momentum is conserved, the contri- 
butions to the spin of the bodies can be computed from 
the orbital contributions: 



Li = -(Gi + G 2 ) 



(6) 



Because we are only interested in the secular evolu- 
tion of the system, we further average the equations of 
motion over the mean anomalies of both orbits (see ap- 
pendix A). The resulting equations are (e.g. Boue and 
Laskar, 2006; Farago and Laskar, 2010): 

Gi = -7(1 - ei)cos/k 2 x ki + 57 (ei • £2)^2 x ei 



> otu cos 6i s,i x ki 



G 2 = -7(1 - el) cos I ki x k 2 + 57 (ei • k 2 )ei x k 2 
- J^ "2i cos £ t §i x k 2 , 
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and 



Lj = -aij cos $i ki x Si - a 2l cos e t k 2 x s, , (10) 

where 
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and 



cos 6i — Sj • ki , cos£i = §i-k2, cos I = ki -k2 ,(14) 

are the direction cosines of the spins and orbits: 9i is the 
obliquity to the orbital plane of the inner orbit, Si is the 
obliquity to the orbital plane of the outer companion, 
and I is the inclination between orbital planes. They 
can also be expressed as 

cos Si = cos I cos Oi + sin I sin Oi cos ipi , (15) 

where ipi is a precession angle. 



compare the theoretical results with the observations, 
as the effect of tides are very small and can only be de- 
tected efficiently after long periods of time. The qual- 
itative conclusions are more or less unaffected, so, for 
simplicity, we adopt here a model with constant Ati, 
which can be made linear (Mignard, 1979; Neron de 
Surgy and Laskar, 1997): 



r[ ~ ri + Ati (wjSj x ri - h) 



(19) 



As for the conservative motion, we can obtain the 
equations of motion directly from equations (4), (5) and 
(6) using Ut instead of U' (see appendix A), that is, 

G 2 = , Gi = -to - U , (20) 



2.2 General relativity correction 

We may add to Newton's equations the dominant con- 
tribution from general relativistic effects. These effects 
are mainly felt by eccentric orbits on close encounters 
between the central bodies and contribute to the gravi- 
tational force with a small correction (e.g. Schutz, 1985) 

2 
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where 1 1 Gi|| = /3iy/ fiiai(l — e\), and c is the speed of 
light. To this order, the dominant relativistic secular 
contribution is on the precession of the pcriapsis, leav- 
ing eccentricity, orientation and semi-major axis of the 
orbit unaffected (Eq. 5): 
3/iini 



ei 



c 2 ai(l - e 2 



■ki x ei . 
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2.3 Tidal effects 

Neglecting the tidal interactions with the external body 
m2, the tidal potential for the inner pair writes (e.g. 
Kaula, 1964): 
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where r^ is the position of the interacting body at a 
time delayed of Ati. The dissipation of the mechani- 
cal energy of tides in the body's interior is responsi- 
ble for this delay between the initial perturbation and 
the maximal deformation. As the rheology of stars and 
planets is badly known, the exact dependence of Ati on 
the tidal frequency is unknown. Many different authors 
have studied the problem and several models have been 
developed so far, from the simplest ones to the more 
complex (for a review see Correia et al, 2003; Efroim- 
sky and Williams, 2009). The huge problem in validat- 
ing one model better than others is the difficulty to 
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The first term in expression (21) corresponds to a 
permanent tidal deformation, while the second term 
corresponds to the dissipative contribution. The pre- 
cession rate of ei about ki is usually much faster than 
the evolution time-scale for the dissipative tidal effects. 
As a consequence, when the eccentricity is constant over 
a precession cycle, we can average expression (23) over 
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the argument of the periapsis and get (Correia, 2009, 
appendix A): 



/ Si + COS 9; ki UJi - \ 

U = -K ini hie,)- — ^ — - / 2 (ei)k! (29) 

\ 2 m l 



3 Secular evolution 

We have presented the equations that rule the tidal evo- 
lution of a hierarchical system of three bodies in terms 
of a vectorial formalism. However, the spin and orbital 
quantities are better represented by the rotation angles 
and elliptical elements. The direction cosines (Eq.14) 
are obtained from the angular momenta vectors, since 
§j = Lj/ ||Lj|| and ki = G»/ ||G»||, as well as the ro- 
tation rate w» = L$ • §j/Cj. The eccentricity and the 
semi-major axis can be obtained from ei = ||ei|| and 
o,i = ||Gi|| 2 /(/3fm(! - ef)), respectively. 



3.1 Conservative motion 

5.1.-/ Cassini states 

The obliquities arc functions of Gi, G2, and L, (Eq. 14). 
Their evolution can be obtained from (Eqs. 7, 8, 10) as 



dcos9i Gi • (§i — cosfljki) L, • (kj 



r// 



|Gi| 
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for the obliquity to the orbital plane of the inner orbit, 
9i. An identical expression could be obtained for the 
obliquity to the outer orbit, £j, replacing Gi by G2. To 
simplify, we may average the equations over the argu- 
ment of the periapsis (appendix A), and the norms of 
||Gi||, HG2II, and ||Lj|| become constant. Thus, 

dcos9 l dki ,. - . dsi ~ 

— 77 — = -7--(s J -cos^ki) + — --(ki-cos^s,) ,(31) 
at at at 
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and 
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at L,; Li 



(33) 



This system has a priori four degrees of freedom asso- 
ciated to the four precession angles (one for each orbit, 
and one per solid body). However, because the total 
angular momentum is conserved, there is only three de- 
grees of freedom. Regular solutions of (Eqs. 31, 33, 33) 



are thus combination of three eigenmodes. More pre- 
cisely, these solutions are composed of a uniform rota- 
tion of all the vectors around the total angular momen- 
tum, and in the rotating frame, vectors describe quasi- 
periodic motions with only 2 proper frequencies (Bouc 
and Laskar, 2006, 2009). The bodies are in an equilib- 
rium configuration, usually called Cassini state, when 
the amplitudes of the quasi-periodic motion vanish. Al- 
though analytical approximations of the solutions can 
be obtained (Boue and Laskar, 2006, 2009), we assume 
that «2i <C a\i <C 7, and also ||Gi|| <C | jGJ-2 1|- Then, 
equation (33) simplifies 

dki 7 /„ . 3 

~~dT " 



1 



cos/k2 x ki 



(34) 



JGiH V ' 2- 

The evolution of Gi is thus independent of Li, and it 
has a uniform precession motion at frequency g around 
the total orbital angular momentum, where 

7 A 3 
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The evolution of the spin-axes is then simpler in the 
rotating frame where ki and k2 are constant. Since 
U\{ 3> a2i, we have for each body 
dsj _ an 

~dt ~ ~Jl7\ 



■ cos 9i ki x §i — g \s.2 x §j , 



(36) 



which leads to 
9i = g sin / sin ipi , 
and 



'■Pi, 



aii 



cos 9i + g cos I ) + g sin I 



r COS (fj 

tan 9j 



(37) 



(38) 



We used the fact that ki and 62 are coplanar. Equation 
(37) shows that the obliquity is oscillating around an 
equilibrium value given by ipi — or n. Stable configu- 
rations for the spin can be found whenever the vectors 
(§i, ki, k2) are coplanar and precess at the same rate 
g (e.g. Colombo, 1966; Peale, 1969). The equilibrium 
obliquities can be found setting (pi — (Eq. 38), which 
provides a single relationship (e.g. Ward and Hamilton, 
2004): 

A, cos 9 t sin 9 t + sin(9 t - I) = , (39) 

where Ai = aij/(||Lj|| g) is a dimensionless parameter. 
The above equation has two or four real roots for 9i, 
which are known by Cassini states. For nearly coplanar 
orbits, we have / ~ 0, and these solutions are approxi- 
mately given by: 



sin/ 



cos/ 



tan Ut±a-J< ±cos {-—)■ (40) 

For a generic value of /, when | Aj | <S^ 1 the first expres- 
sion gives the only two real roots of equation (39). On 
the other hand, when | Aj | 3> 1 , we have four real roots 
approximately given by expressions (40). 
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3.1.2 Lidov-Kozai cycles 

The variations in the mutual inclination between the 
inner orbit and the orbit of the external companion 
can be obtained from the direction cosine (Eq.14) in a 
similar way as the obliquity (Eq. 30): 



dcosl Gi • (k2 — cos/ki) G2 ■ (ki — cos Ik^) 

dt fiGJ + \^\\ 

Since ||Gi|| <C || G2 1| , we obtain from expression (7) 
dcosl 5 ~/e\ 



(41) 
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2||Gi 
^-^ an cos 9. 



^ IIGi 



cos / sin / sin 1w\ 
(§i x ki) k 2 



(42) 



where w\ is the argument of the periapsis of the inner 
orbit, that is, the angle between the line of nodes of the 
two orbits and the periapsis of the inner orbit. 

On the other hand, the variations in the eccentricity 
are easily obtained from the Laplace- Runge-Lenz vector 
(Eq.9): 

gi -ei _ 5 7(1 - ef)ei ._ 
~ 2" 



ei 



■ sin / sin 2u?i 



(43) 



ei -i ||Gi| 

Thus, combining expressions (43) and (43) and ne- 
glecting the contributions from the rotational flattening 
(terms in aij), we get 

rfcosJ - e 2h cosI (44) 

dt "(l-ef) COS7 ' l44j 

which can be integrated to give 



1 — e\ cos I = hi = Cte 



(45) 



where hi is constant in absence of tides. The above re- 
lation shows that an increase in the eccentricity of the 
inner orbit must be accompanied by a decrease in the 
mutual inclination and vice-versa. It corresponds to a 
happy coincidence, as called by Lidov and Ziglin (1976), 
which results from the fact that in the quadrupolar ap- 
proximation the potential energy (Eq. 2) does not de- 
pend on argument of the periapsis of the external body 
za2 (e.g. Farago and Laskar, 2010). As a consequence, 



The variations in the inclination and eccentricity 
(Eqs. 43, 43) become zero when sin2-n7i = 0, that is, 
for cos2zui = ±1. When cos2o7i = —1 (vj\ = ±7r/2) it 
is possible to show that w = has a solution if hi < 
a/375 (Lidov, 1962; Kozai, 1962). Thus, when the mu- 
tual inclination is greater than arccos \/3/5 ~ 39.23° 
we have a libration regime for the periapsis about w\ = 
±7r/2. In this regime, one can observe large variations 
in both I and ei, known by Lidov-Kozai cycles. If the 
inner orbit is initially circular, the maximum eccentric- 
ity achieved is given by e\ = y/l — (5/3) cos 2 /. 

Lidov-Kozai cycles persist as long as the perturba- 
tion from the outer body is the dominant cause of pre- 
cession in the inner orbit. However, additional sources 
of precession, such as general relativity or tides, can 
compensate the libration mechanism and suppress the 
large eccentricity/inclination oscillations (e.g. Migaszewski 
and Gozdziewski, 2009). For hi > \/3/5 the periapsis 
of the inner orbit is always in a circularization regime, 
so there is only small variations in the eccentricity and 
inclination. 



3.2 Tidal evolution 

3.2.1 Spin evolution 

The variation in the body's rotation rate can be com- 
puted from equation (29) as uii = Lj • Sj/Cj, giving 
(Correia and Laskar, 2010a): 

Km (,. l + cos^Wi \ 

Wi= 77- /i(ei) /2(ei)cos0i .(47) 

For a given obliquity and eccentricity, the equilib- 
rium rotation rate, obtained when w, = 0, is attained 
for: 



the conjugate Delaunay variable || G2 1| = /?2\/M2 a 2(l — e|) 
is constant, and, as the semi-major axis are constant in 
the secular problem, e2 is also constant. Indeed, when 
we consider only the orbital contributions, the total an- 
gular momentum is 
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2||Gi|| ||G 2 ||cosJ = Cie 



(46) 



Since || G-2 1| is constant and ||Gi|| <C HG2II, it remains 
||Gi|| cos/ = Cte, and we retrieve the result given by 
equation (45). This is no longer true for the octopole 
or higher order approximations (e.g. Laskar and Boue, 
2010). 



/ 2 (ei) 2cos( 



(48) 



m /i(ei) 1 + cos 2 6»i 

Notice, however, that the above expression is only valid 
for constant eccentricity (and also constant semi-major 
axis) , or at least if tidal effects modify the eccentricity 
faster than other effects. Indeed, when the eccentric- 
ity is forced by the gravitational perturbations from a 
companion body, the limit solution of expression (47) 
is no longer given by equation (48), but more generally 
(Correia and Laskar, 2004, 2009): 

t 
f 2 (e 1 (T)) cos 6i(T) 9i (T) dr , (49) 
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It corresponds to a solution that pursues the instan- 
taneous equilibrium rotation for a given eccentricity 
(Eq.48), but delayed and with a smaller amplitude, de- 
pending on the relative strength of tidal effects. The 
stronger these effects are, the shorter is the time delay 
and closer are the amplitudes to expression (48). 

In turn, the dissipative obliquity variations are com- 
puted by substituting equation (29) in (30) with ||Lj|| <C 
||Gi||, giving: 

'' h ' ni sm$ i (f 1 (e 1 )cose i ^-f2(ei)) . (51) 



CiUJi \ J *~ y ''' ' 2ni 

Because of the factor ni/wj in the magnitude of the 
obliquity variations, for an initial fast rotating body, 
the time-scale for the obliquity evolution is longer than 
the time-scale for the rotation rate evolution (Eq.47). 
As a consequence, it is expected that the rotation rate 
approaches its equilibrium value (Eq.48) earlier than 
the obliquity. Replacing equation (48) in (51), we have 
for constant eccentricity: 



K l n 1 
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h(ei) 



sin 9j 



1 + cos 2 6i 



(52) 



We then conclude that, although the initial behavior of 
the obliquity depends on the initial rotation rate, tidal 
effects always end by decreasing the obliquity, since Qi < 
0. Thus, the final obliquity tends to be captured in a 
small obliquity Cassini state (Eq.40), that is, 

sin/ CiUJij ( _ 3 



Xi 
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(53) 



3.3 Orbital evolution 

The variations in the norm of the orbital angular mo- 
mentum can be computed directly from expression (29), 
since Gi = — Lo — La: 

| ||Gl|| = -£>•]<! 
i 

= ^tfini(A(e 1 )cos0^-/ 2 (ei)) . (54) 

The variations in the eccentricity are easily obtained 
from the Laplace- Runge-Lenz vector (Eq.21): 

ei -ei 
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ei 



= £ Iff (i^ ei ) cos ^ ^) ei ' ^ 

while the semi-major axis variations are obtained from 
the eccentricity and the norm of the orbital angular 
momentum: 

hi 2eiei 2 Gi • Gi 

~i~(l-ej) + \\Gif 
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(56) 



The inner orbit can either expand or contract, de- 
pending on the initial spin of the bodies. Considering 
for simplicity dissipation in only one component, for 
fast initial rotating rates (ui 3> ni), the semi- major 
axis and the eccentricity usually increase, except for 
retrograde spins ($i > 7t/2). The ratio between or- 
bital and spin evolution time-scales is roughly given 
by Ci/(mia\) <C 1, meaning that the spin achieves an 
equilibrium position much faster than the orbit. Thus, 
as the rotation rate decreases, the increasing tendency 
in the orbital parameters is reversed when d ||Gi || jdt — 
(Eq.54), 

— cos Qi = , (57 

n\ fi(ei) 

that is, when the rotation rate is close to its equilib- 
rium value (Eq. 48). After an initial increase in the semi- 
major axis and in the eccentricity, we can then always 
expect a contraction of the inner orbit until it becomes 
completely circularized. The final evolution of the sys- 
tem is then achieved when ei = 0, uji = rii (Eq.48) and 
$i ~ -sin //A* (Eq.53). 

The variations in the mutual inclination can be ob- 
tained using expression (41) with Gi = — J^i ^t (Eq- 29). 
We have then: 
d cos / ^— v K iUi 



dl 



^2 lid 



■ f i(ei) (cos Ei — cos I cos 8i) , (58) 



or, making use of expression (15), 
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(59) 



Just after the spin of one body is trapped in a small 
obliquity Cassini state, ifi — and sin^ is constant 
and given by expression (53). Thus, 

dl s—^ Kioji 



dt ^ 2Ai IIGi 



r/i(ei)sin/ 



(60) 



4 Application to exoplanets 

In this section we apply the model described in Sect. 2 
to three distinct situations of exoplanetary systems: 
HD 80606 (inner restricted problem), HD 98800 (outer 
restricted problem), and HD 11964 (intermediate non- 
restricted problem). We numerically integrate the set 
of equations (7), (8), (9), (10) and (17) for the conser- 
vative motion, together with equations (20), (21) and 
(23) for tidal effects. 

There exist many systems containing a "hot Jupiter" 
in a wide binary for which one could apply the present 
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model to illustrate the tidal migration combined with 
Lidov-Kozai cycles. However, we prefer to reproduce the 
results for HD 80606 in order to compare our results 
(obtained with a non-restricted model) with previous 
studies by Wu and Murray (2003) and Fabrycky and 
Tremaine (2007) (obtained with a restricted model). 

The stability of a disc around the multi-binary HD 98800 
system has been recently studied by Verrier and Evans 
(2009) and Farago and Laskar (2010), which have iso- 
lated two different regimes for the trajectories of the 
disc. We then apply our model to an hypothetical planet 
around the binary stars and show how transitions be- 
tween the two regimes are possible, and how particles 
in initial prograde orbits may end in retrograde orbits. 

We finally use our model to study the HD 11964 
system, which is a hierarchical planetary system com- 
posed of two planets, the inner one in the Neptune-mass 
regime and the outer one similar to Jupiter. There is 
no determination of the mutual inclination between the 
two planets, but the system is stable for very large val- 
ues, so we analyze its behavior for different situations. 



Table 1 Initial parameters for the HD 80606 system (Naef 
et al, 2001; Eggenberger et al, 2004; Wu and Murray, 2003). 



4.1 HD 80606 

In current theories of planetary formation, the region 
within 0.1 AU of a protostar is too hot and rarefied 
for a Jupiter-mass planet to form, so "hot Jupiters" 
likely form further away and then migrate inward. A 
significant fraction of "hot Jupiters" has been found in 
systems of binary stars (e.g. Eggenberger et al, 2004), 
suggesting that the stellar companion may play an im- 
portant role in the shrinkage of the planetary orbits. 
In addition, close binary stars (separation compara- 
ble to the stellar radius) are also often accompanied 
by a third star. For instance, Tokovinin et al (2006) 
found that 96% of a sample of spectroscopic binaries 
with periods less than 3 days has a tertiary compo- 
nent. Indeed, in some circumstances the distant com- 
panion induces tidal interactions in the inner binary 
by means of the Lidov-Kozai mechanism (Sect. 3.1.2), 
causing the binary semi-major axis to shrink to the 
currently observed values (e.g. Eggleton and Kiseleva- 
Eggleton, 2001). The same mechanism has been subse- 
quently proposed to be at the origin of "hot Jupiters" 
as an alternative to migration in a disk (e.g. Wu and 
Murray, 2003). 

The HD 80606 system is composed of two Sun- like 
stars in a very wide orbit (02 ~ 1000 AU) (Eggenberger 
et al, 2004), and a short-period planet in a very eccen- 
tric orbit (ai = 0.45 AU and e x = 0.92) (Naef et al, 
2001). Because at periapsis the distance to the main 
star is only 0.036 AU, the orbit of the planet is still un- 
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dergoing tidal evolution, and is thus a perfect example 
to test our model. 

In order to compare our results with the previous 
studies of Wu and Murray (2003) and Fabrycky and 
Tremaine (2007) we use the same initial conditions for 
the planet and the stellar companions: the planet is 
initially set in a Jupiter-like orbit with a\ = 5AU, 
ei = 0.1, and / = 85.6°, while the stellar companion 
is supposed to be a Sun-like star at a2 = 1000 AU, 
and e2 = 0.5 (Table 1). In Figure 1 we plot an ex- 
ample of combined tidal-Kozai migration of the planet 
HD 80606 b. 

Prominent eccentricity oscillations are seen from the 
very beginning and the energy in the planet's spin is 
transferred to the orbit increasing the semi-major axis 
for the first 0.1 Gyr (Eq. 47). As the equilibrium ro- 
tation is approached around 0.3 Gyr (Eq. 49) the tidal 
evolution is essentially controlled by equations (55) and 
(56), whose contributions are enhanced when the eccen- 
tricity reaches high values. The semi-major axis evo- 
lution is executed by apparent "discontinuous" transi- 
tions precisely because the tidal dissipation is only effi- 
cient during periods of high eccentricity. As dissipation 
reduces the semi-major axis, periapsis precession be- 
comes gradually dominated by general relativity rather 
than by the third body and the periapsis starts circu- 
lating as the eccentricity approaches to near 0.5 Gyr. 
Tidal evolution stops when the orbit is completely cir- 
cularized. The final semi-major axis is estimated to 
about a,f = ai(l — e\) ~ 0.07 AU, which corresponds 
to a regular "hot Jupiter" (Correia and Laskar, 2010b). 

The angle between the spin axis of the planet and 
its orbit (denoted by #1) is quickly brought by tides 
to a small obliquity Cassini state (Eq.53). On the other 
hand, the angle between the spin axis of the star and the 
orbit of the planet (denoted by 8 n ) is not tidally evolved, 
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Fig. 1 Possible evolution of the planet HD 80606 b initially in an orbit with a\ — 5 AU, e± = 0.1, and / = 85.6° (Table 1). We 
show the evolution of the eccentricity e\ (a), semi-major axis a± , and periapsis a± (1 — ei) (b), mutual inclination /, obliquity of 
the star f?o, and obliquity of the planet 9± (c), and orbital period P or b, rotation period of the star P r ot,0j and rotation period 
of the planet P ro t,i (d). Results plot here are identical to those in Wu and Murray (2003) and Fabrycky and Tremaine (2007), 
obtained with a restricted model. 



since the dissipation within the star is much smaller 
than the dissipation within the planet (fc2 Ato <C fei At\ , 
Tablet). As a consequence, capture of the spin of the 
star in a Cassini state does not occur during the time 
length of the simulations, but one can observe large os- 
cillations corresponding to the variations in the orienta- 
tion of the orbital plane of the planet. Indeed, as long as 
«io "C 7, the angle between the spin of the star and the 
outer companion (denoted by Eq) is approximately con- 
stant (Eq.36), and thus |e - I\ < 6 < £q + I (Eq. 15). 
In our simulation e = 85.6° - 10° = 75.6° (Table 1), so 



we approximately observe that 10 c 
urelc). 



< 6» < 161° (Fig- 



As the semi-major axis decreases, the precession 
of the planet's orbit becomes progressively dominated 
by the equatorial bulge of the star, so that aio ~ 7 
around 2.8 Gyr (Eq. 33). From that point, the orbit of 
the planet precesses around the spin of the star and 6$ 
becomes constant, retaining a value between 10° and 
160°. One then expects to observe a misalignment be- 
tween the spin of the star and the orbit of the planet. In 
the simulation shown in Figure 1, we obtain 9q « 60°, 
but this value depends on the initial configuration of 
the system. 

In order to get a statistical distribution of the fi- 
nal misalignment we have integrated a series of 40 000 
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Table 2 Initial parameters for the HD 98800 B system (Tor- 
res et al, 1995; Tokovinin, 1999). 



Fig. 2 Histogram of the final distribution of the misalign- 
ment angle 9q ■ We have integrated a series of 40 000 systems 
with the same initial conditions from Table 1 except for the 
obliquity 0q = 0°, and inclination /, which range from ±84.3° 
to 90° . We observe two pronounced peaks of higher probabil- 
ity around 6q r; 53° and Gq r; 109°, which is consistent with 
the observations of the Rossiter-McLaughlin anomaly for the 
HD 80606 system (Pont et al, 2009). 
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systems with the same initial conditions from Table 1 
except for the obliquity of the star 9 = 0° , and for the 
mutual inclination /, which ranged from ±84.3° to 90° 
(on an evenly spaced grid of —0.1 < cos/ < +0.1). A 
histogram with the results of these simulations is given 
in Figure 2, with a bin width of 1°. Our results do not 
differ much from the histogram obtained by Fabrycky 
and Tremaine (2007), except that we observe two pro- 
nounced peaks of higher probability around 6q f» 53° 
and 6*o ~ 109°. They performed 1000 simulations dis- 
tributed in bins widths of 10°, which partially explains 
this difference, but the main reason is the fact that they 
adopted for the rotation of the star P ro t,o = 10 h, while 
we used a value closer to the Sun's rotation period. Be- 
cause we executed so much simulations our histogram 
is close to the probability density function distribution 
for the misalignment. Observations of the Rossiter- 
McLaughlin anomaly for the HD 80606 star reinforced 
the hypothesis of spin-orbit misalignment in this sys- 
tem (alignment excluded at > 95% level), with a pos- 
itive median projected angle of 50° (Pont et al, 2009). 
These observations are in perfect agreement with our 
simulations. 

Although we used the full quadrupolar problem, 
while previous studies ( Wu and Murray, 2003; Fabrycky 
and Tremaine, 2007) used the inner restricted problem 
(where the outer orbit is considered constant), we re- 
trieve similar results and then verify that their approx- 
imation is appropriate for this specific situation. 



4.2 HD 98800 

HD 98800 is an interesting and unusual system of four 
lOMyr old post T-Tauri K stars: two spectroscopic bi- 
naries A and B in orbit about one another (Torres et al, 



1995; Kastner et al, 1997). Moreover, it has a large in- 
frared excess attributed to a circumbinary disc around 
the B pair (Koerner et al, 2000; Furlan et al, 2007). 
Substantial extinction towards this pair suggests that 
it is observed through some of this material (Tokovinin, 
1999). An apparent absence of CO molecular gas in the 
disc combined with infrared spectrum modeling indi- 
cate that this is a T-Tauri transition disc that is just 
reaching the debris disc stage, with a collisional cascade 
having been recently initiated (e.g. Furlan et al, 2007). 
The orbits of the stars are all highly eccentric and in- 
clined, creating a dynamical environment unlike almost 
all other known debris discs (Tokovinin, 1999). The dust 
disc is generally agreed to be an annulus around the 
B binary, but the exact structure varies from model to 
model. Koerner et al (2000) estimate a coplanar nar- 
row ring outwards of about ~ 5 AU from the two stars. 
However, Boden et al (2005) argue that the line of sight 
extinction means that the disc cannot be coplanar un- 
less it is very flared. 

Verrier and Evans (2009) investigated numerically 
the stability of a family of particles at large inclinations 
around the B binary, which remain stable even under 
the perturbation of an outer third stellar companion. 
The results show that the Lidov-Kozai mechanism of 
the outer star is disrupted by a nodal libration induced 
by the inner binary pair on a shorter time-scale. An 
analytical study by Farago and Laskar (2010) confirmed 
that equilibria at large inclination circumbinary orbits 
exist for the outer restricted quadrupolar problem. This 
raises the possibility that planets and asteroids with 
large inclination may survive in multi-stellar systems. 

In this section we test the evolution of a Jupiter- 
like mass planet around the B binary in an orbit that 
is presently occupied by the dust disc. Since we aver- 
age the orbit of the planet over the mean longitude, the 
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Fig. 3 Energy levels in the (7 sin to, I cos to) plane for different values of the eccentricity of the B binary HD 98800, e\ = 0.5 
(a), ei = 0.785 (b), and ei = 0.9 (c) (Farago and Laskar, 2010). The dot marks the initial position of the outer orbit (Table 2). 



orbital evolution of this planet can also be regarded as 
the orbital evolution of the particles in the disc. The 
HD 98800 system is still very young, and the B binary 
stars are very close to each other, so one can expect 
that the two components will undergo significant tidal 
evolution throughout its life. In order to speed-up the 
simulations we assume that both stars experience in- 
tense tidal dissipation, with At = 10 2 s (Table 2). Be- 
cause the gravitational effects of the A binary do not 
destabilize the system, we do not consider its presence 
in the simulations. 

The inclination of the B binary with respect to the 
line of sight is estimated to be about 23° (Boden et al, 
2005). The inclination of the disc with respect to the 
B binary is unknown, but one can contest its copla- 
narity, since substantial extinction of the light implies 
that the system is observed through some of this ma- 
terial (Tokovinin, 1999). Therefore, we assume that the 
inclination is about 20° (Table 2). The argument of the 
periapsis and the longitude of the node were also de- 
termined with respect to the plane of the sky (Boden 
et al, 2005), but again, we miss the value of these two 
quantities measured with respect to the plane of the 
disc. As a consequence, w\ is a free parameter for the 
planet/disc in our simulation. The choice of this pa- 
rameter is critical, as it places the initial system in a 
different regime of libration, and it also determines the 
libration amplitude (Fig. 3b). Indeed, we can differen- 
tiate two kinds of regimes (Farago and Laskar, 2010): 
closed trajectories where the orbital angular moment 
of the planet k2 precesses around the orbital angular 
momentum of the binary ki (or its opposite — ki); or 
closed trajectories where the orbital angular momen- 
tum of the planet precesses around the direction of the 
periapsis of the binary ei (or its opposite — ei). In the 
first situation the inclination is strictly inferior to 90° 
(or strictly superior to 90°, corresponding to a retro- 



grade orbit), while in the second case the inclination 
oscillates around ±90°. 

In our simulation, we use w\ = 82°, as this value 
places the initial orbit of the planet at the edge between 
the two different regimes of behavior (Fig. 3b) . Initially, 
k2 is precessing around ki and the mutual inclination 
oscillates between 20° and nearly 90°. As the system 
evolves, the norm of the orbital angular momentum of 
the planet remains constant, but the norm of the or- 
bital angular momentum of the binary is modified by 
exchanges with the rotational angular momentum of 
the stars. Since we have considered wi 3> w (Table 2) 
in our example the initial exchanges are mainly between 
the binary's orbit and the spin of the star 1. 

During the first stages of the evolution wi ^$> n\ and 
the norm of the orbital angular momentum of the bi- 
nary increases (Eq. 54), as well as the eccentricity and 
the semi-major axis (Eqs. 55, 56). As a consequence, the 
separatrix between the two regimes contracts (Fig. 3c). 
Since the norm of the orbital angular momentum of 
the planet remains constant (the planet is too far to 
undergo dissipation), the orbit of the planet crosses 
the separatrix and switches to the regime of preces- 
sion around the direction of the periapsis of the binary 
ei. In our simulation this transition occurs shortly after 
80 Myr since we placed the initial system very close to 
the separatrix (Fig. 4a). 

According to expression (51) the obliquity of the 
star 9i also increases, but the rotation rate decreases 
(Eq. 47). Around 2 Gyr, the obliquity begins to decrease 
and the orbit of the binary initiates its contraction. In 
general, this contraction is much slower than the pre- 
vious expansion because the spin is already near an 
equilibrium position and the major source of dissipa- 
tion becomes the orbital energy. As the eccentricity is 
damped, the separatrix between the two regimes ex- 
pands and catches again the planet's orbit (Fig. 3b). 
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Fig. 4 Possible evolutions for a rri2 ~ 1O -3 M0 planet or disc at a-2 = 5.2 AU and e-2 = 0.5, initially in a prograde orbit around 
the binary HD 98800 B (Table 2). During the first evolutionary stages the orbit of the binary expands and the planet changes 
of regime after about 80 Myr (a). As the binary orbit shrinks, the planet changes of regime again (somewhere between 25 
and 30Gyr), but depending when the planet crosses the separatrix between regimes, its orbit may become retrograde (b) or 
prograde (c). 



After the second separatrix transition, the orbital an- 
gular momentum of the planet precesses again around 
the orbital angular momentum of the binary or its op- 
posite ±ki , depending on the w\ value at the time the 
transition occurs. We have performed two different sim- 
ulations, with a 10~ 2 difference in the initial obliquity 
of the star 9q, and each one led to the a different final 
evolution scenario (Fig. 4b, c). 

Subsequent tidal evolution of the binary orbit tends 
to circularize it. When the eccentricity reaches zero, the 
oscillations in the mutual inclination are also damped 
and the two angular momenta precess at constant incli- 
nation and speed (Fig. 3a). Tidal evolution in the inner 
binary is thus a very efficient mechanism of transform- 
ing initial prograde (or retrograde) orbits in retrograde 
(or prograde) orbits, although the norm of the angular 
momentum of the planet does not change. 



4.3 HDfl964 

The planetary system around HD 11964 is composed of 
two planets with minimum masses mi = 25 M ffi (planet 



c) and m 2 = 0.62 M Jup (planet 6), at a\ = 0.229 AU and 
a 2 = 3.16AU, respectively (Butler et al, 2006; Wright 
et al, 2009). Veras and Ford (2010) performed extensive 
n-body simulations for this system and concluded that 
it is always stable for mutual inclinations / < 60° (or 
/ > 120°), stable up to 85% for / < 75° (or / > 105°), 
and stable up to 25% for 75° < / < 105°. It is a system 
for which stability is possible for a wide range of mu- 
tual inclinations and thus perfect to apply our secular 
model. 

Since the ratio between the semi-major axis of the 
two planets is aij a\ ~ 14, the quadrupolar approxima- 
tion is well suited, although octopole order terms in the 
development of the potential energy (Eq. 2) may cause 
variations and exchanges between the eccentricities of 
both planets. In order to test the quality of our secu- 
lar model in the case of hierarchical planetary systems 
of this kind, we performed some direct n-body numeri- 
cal simulations for the HD 11964 system and compared 
with the quadrupolar secular model. Results for the ini- 
tial parameters in Table 3 are shown in Figure 5. We 
observe slightly additional variations in the amplitudes 
and precession rates of the eccentricity and inclination, 
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Table 3 Initial parameters for the HD 11964 system (Butler 
et al, 2006; Wright et al, 2009). 
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but the general long-term behavior of the system re- 
mains essentially the same. 

Though the present age of the star is estimated to 
be about lOGyr (Saffe et al, 2005), we start the sys- 
tem as if it has been recently formed. We assume the 
present orbits as the original ones, and an initial rota- 
tion period of the inner planet of about 10 h (Table 3). 
Because tidal effects on the inner planet are strong, the 
choice of the initial spin is irrelevant, as it evolves to 
an equilibrium configuration in about one million years 
(Figure 6a). Nevertheless, this procedure allows us to 
understand simultaneously the initial behavior of the 
system and also its future evolution. In addition, we 
arbitrarily assume a mutual inclination / = 30°. This 
inclination value sets the system well outside a coplanar 
configuration, but yet below the critical value of about 
39° that allows Lidov-Kozai cycles. 

During the first million years (Figure 6a), the ro- 
tation rate of the inner planet and its obliquity are 
brought to their equilibrium positions, where one be- 
lieve that the spin of the planet can be found today. The 
rotation rate is slowed down until it reaches an equilib- 
rium near a>i/ni « /2(ei)//i(ei), that is, P ro t,i ~ 25d 
(Eq. 49) , while the obliquity, after an initial increase 
(since the initial rotation rate is fast (Eq.51)), quickly 
drops until it is captured in the small obliquity Cassini 
state 6»i w sin//Ai « 0.018° (Eq.53). Because the ec- 
centricity and the inclination are oscillating, so does the 
equilibrium rotation rate (Eq. 48) and the equilibrium 
obliquity (Eq.53). 

During the first two million years (Figure 6a), the 
eccentricity and the inclination of the inner orbit do not 
experience any substantial secular modification (apart 
the periodical ones) , since tidal effects are much less ef- 
ficient over the orbit than over the spin. However, when 
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Fig. 5 Evolution of the HD 11964 c inclination (top) and ec- 
centricity (bottom) during 1 Myr, starting with the orbital 
solution from Wright et al (2009) and / = 30° (Table 3). The 
solid lines curves are the values obtained with the quadrupo- 
lar secular model, while the dashed green lines are the com- 
plete solutions obtained with a n-body numerical model. 



we follow the system over Gyr time-scales one can ob- 
serve those variations occurring (Figure 6b). The most 
striking effect is an initial increase in the eccentricity 
of the inner planet, that is accompanied by a signifi- 
cant reduction in the mutual inclination between the 
orbits of the two planets. This behavior results from a 
combination of tidal effects, rotational flattening, and 
the gravitational perturbations from the outer planet, 
which tend to conserve the quantity h\ = y/l 



e\ cos / 



(Eq. 45). Tidal effects alone also decrease the eccentric- 
ity (Eq. 55), so when its contribution becomes dominat- 
ing the eccentricity is progressively damped to zero. 

The exchanges between eccentricity and inclination 
(Fig. 5) are more significant for large values of the ini- 
tial mutual inclination. The maximum value for the ec- 
centricity oscillations is then higher and therefore the 
planet experiences stronger tidal effects. As a result, 
the orbit of the inner planet evolves in a shorter time- 
scale. In Figure 7 we plot the simultaneous evolution 
of the inclination and eccentricity of the inner orbit 
for different initial mutual inclinations, I. For / < 30° 
there is almost no reduction in the inclination, and the 
eccentricity damping is very slow (Figure 7a). On the 
other hand, for / > 40° we observe a reduction of more 
than 20° in the inclination, which is accompanied by an 
initial increase in the eccentricity, followed by a rapid 
damping to zero (Figure 7c, d). 

For large values of the initial inclination, the system 
evolves much faster to its final configuration, so the or- 
bit becomes circular in a shorter time-scale. Since the 
present orbit of the inner planet still presents an ec- 
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Fig. 6 Long-term evolution of the planet HD 11964 c for an initial mutual inclination of 30° (Table 3). During the first million 
years, the rotation rate is slowed down into its equilibrium position (Eq. 48) and the obliquity is trapped in a small obliquity 
Cassini state (Eq. 53). As the system evolves, tidal effects decrease the mutual inclination of the planetary orbits, while the 
eccentricity increases until around 2 Gyr, time after which it is damped to zero. The rotation rate closely follows the eccentricity, 
and the obliquity also decreases, as the equilibrium in a Cassini state depends on the inclination of the perturber. 
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Fig. 7 Long-term evolution of the inclination (top) and eccentricity (bottom) of the planet HD 11964 c for different values of 
the initial mutual inclination: 15° (a), 30° (b), 40° (c), and 60° (d) (Table 3). As the initial inclination increases, so does the 
amplitude of the eccentricity oscillations, resulting that tidal effects are enhanced, and the evolution time-scale shorter. 
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Fig. 8 Time needed to circularize the orbit of the planet HD 11964 c (ei < 0.01) as a function of the initial mutual inclination 
and for different values of the initial eccentricity (ei = 0.1, 0.2, 0.3, and 0.4). The system is 10 Gyr old (Saffe et al, 2005), so 
we conclude that initial inclinations between 40° and 140° were not possible, as the present eccentricity is still 0.3. 



centricity close to 0.3, one may assume that the initial 
inclination could not have been excessively large. In or- 
der to estimate an upper limit for the initial inclination, 
we can run some simulations and check which of those 
can evolve to the present configuration within 10 Gyr, 
the estimated age for the system. 

However, the choice of the initial orbit is not easy, 
since many different initial configurations can bring the 
planet to the present orbit, depending on the tidal damp- 
ing coefficients (fc2 and At) and also on the initial eccen- 
tricity, inclination, and semi-major axis. It is out of the 
scope of this study to explore exhaustively the initial 
configurations of the HD 11964 system, but rather to il- 
lustrate some interesting evolutionary scenarios. As we 
just have seen, the present eccentricity combined with 
a large inclination can later evolve to an identical ec- 
centricity but with smaller inclination. The semi-major 
axis will be slightly smaller than today, but this pa- 
rameter essentially plays over the evolution time-scale. 
Thus, the present system can be a representation of the 
initial system and we adopt it for simplicity. In Fig- 
ures we used the initial conditions listed in Table 3, 
but modified the initial inclinations and eccentricities. 
We observe that inclinations in the interval [40°, 140°] 
circularize the orbit (ei < 0.01) in less than 10 Gyr, so 
they can be discarded. 



5 Conclusion 

Many multi-planet systems have been reported in hier- 
archical configurations. For the most part, their mutual 
inclinations are unknown, but the fact that they exhibit 
significant values in the eccentricities led to think that 
the inclinations can also be large. In addition, very of- 
ten the innermost planet in these systems is very close 



to the star and undergoes tidal evolution. Here we have 
studied this particular sub-group of multi-planet sys- 
tems. Using a very general and simplified averaged vec- 
torial formalism, we have shown that inclined hierar- 
chical planetary systems undergoing tidal dissipation 
can evolve in many different and sometimes unexpected 
ways. 

We re-analyzed the case of HD 80606, a system where 
the planet migrated inwards by a combined tidal-Kozai 
mechanism, confirming previous results. We looked at 
the behavior of a planet (or disc) around the binary 
stars HD 98800 B and showed that initial prograde or- 
bits may become retrograde and vice-versa, only be- 
cause of tidal migration within the binary stars. Fi- 
nally, we studied the regular 2-planet HD 11964 system 
and showed that tidal dissipation combined with grav- 
itational perturbations may lead to a decrease in the 
mutual inclination, and a fast circularization of the in- 
ner orbit. 

We have chosen the above three examples, as they 
are representative of the diversity of behaviors among 
inclined hierarchical systems. Many other systems are 
awaiting to be studied. The fact that we use average 
equations for both tidal and gravitational effects, makes 
our method suitable to be applied in long-term studies. 
It allows to run many simulations for different initial 
conditions in order to explore the entire phase space 
and evolutionary scenarios. In particular, it can be very 
useful to put constraints on the inclinations and dissipa- 
tion ratios of hierarchical planetary systems. Our study 
can also be extended to systems of binary stars, and to 
planet-satellite systems. 
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A Averaged quantities 

For completeness, we gather here the average formulae that 
are used in the computation of secular equations. Let F(r,r) 
be a function of a position vector r and velocity r. Its averaged 
expression over the mean anomaly (M) is given by 

1 f 27V 
(F)m=7T / f ( r » dM ' ( 61 ) 

2n Jo 

Depending on the case, this integral is computing using the 
eccentric anomaly (E), or the true anomaly (v) as an inter- 
mediate variable. The basic formulae are 



dM = -&E = 



=dv 



a 2 V^ 

r = a(cos E — e) e + a\J\ — e 2 (sin Ejkxe, 
r = r cos v e + r sin »kxe, 

: k X (r + e) , 



A 



r = a(l — e cos B) 



a(l-e 2 ) 
1 + e cos v 



(62) 



where k is the unit vector of the orbital angular momentum, 
and e the Laplace- Runge-Lenz vector (Eq. 3). We have then 



l\ = I 

r 3 / a 3 (l-e 2 ) 3 / 2 

and 



- 1 (i _ k * k ^ 

r 5 J 2a 3 (l-e 2 ) 3 / 2 V ' ' 



which leads to 
-P 2 (f-u) 



2a 3 (l-e 2 ) 3 / 2 
for any unit vector u. In the same way 

(, 2 ) = „ 2 (l + V) , 



P 2 (k-u) , 



and 



.1- 



( r <r)=„ 2 (l-k<k)+±a 2 



(63) 



(64) 



(65) 



(66) 



(67) 



give 

(r 2 P 2 (r-u)) 



(l-e 2 )P 2 (k-u)-5e 2 P 2 (e-u) . (68) 



The other useful formulae are 



le ) = — M e ) < 



8 / a8vT-e 2 



/2(e) 



(69) 

(70) 



r r 



2a D 
i' \ 5 1 



V '- e!, /,(e)(l-k'k) + < , H/ ; + l q/9 e'e,(71) 



2.7./, e 2 



' a 7 V^- 



4a 6 (l-e 2 ) 9 / 2 
/4(e)e , (72) 



7 1 



r 10 / 2a 9 (l-e 2 
(r • r)r 



■/s(e)€ 



2a 7 yr 



= / 5 (e)k X e , 



(73) 



(74) 



where the /i(e) functions are given by expressions (24) to 
(28). 

Finally, for the average over the argument of the periapsis 
(va), we can proceed in an identical manner: 

<e 4 e) ro = ^ r e *e^=^(l-k'k) , (75) 



which gives 

<(e-u)e> w = y(u-(k.u)k). (76) 
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